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Abstract We present formulae for accurate numerical conversion between func- 
tions represented by multiwavelets and their multipole/local expansions with re- 
spect to the kernel of the form, e~ Ar /r (cf. [7 ). The conversion is essential for the 
application of fast multipole methods for functions represented by multiwavelets. 
The corresponding separated kernels exhibit near-singular behaviors at large A. 
Moreover, a multiwavelet basis function oscillates more wildly as its degree in- 
creases. These characteristics in combination render any brute-force approach 
based on numerical quadratures impractical. Our approach utilizes the series ex- 
pansions of the modified spherical Bessel functions and the Cartesian expansions 
of solid harmonics so that the multipole-multiwavelet conversion matrix can be 
evaluated like a special function. The result is a quadrature-free, fast, reliable, 
and machine precision accurate scheme to compute the conversion matrix with 
predictable sparsity patterns. 

Keywords Multipole expansion ■ Multiwavelet • Screened Coulomb Potential 
1 Introduction 

Multiwavelets [IJ[2] developed originally by Alpert generalize Haar wavelets with 
piecewise polynomial scale functions up to any given degree, hence, enjoy higher 
order of accuracy in the representation of sufficiently smooth functions. A de- 
tailed discussion on the sparse representation of differential operators and expo- 
nential operators for evolution equations can be found in [2]. Recent advances in 
multi-resolution algorithms for integral operators can be found in various articles 
including 13,4,5,6 , which are based on the dimensional kernel separation via rep- 
resentation of kernels by weighted sum of Gaussians, thus, applicable to a variety 
of kernels in arbitrary dimensions. 
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In this paper, we focus ourselves on a more traditional, but very popular and 
well-studied fast convolution algorithm - the fast multipole method. Our goal 
is to establish the connection between multiwavelet representation and the fast 
multipole method. The problem can be reduced to finding the multipole expansion 
of multiwavelet basis functions. 

Following one of the most recent version of the fast multipole algorithm pre- 
sented in 0, we begin with the kernel, 

P -H«>-y\\ 

where A is a non-negative real number. This kernel is the fundamental solution of 
the linear differential operator, 

A — A 2 , (2) 

which appears in various applications involving damped Coulomb forces. The re- 
sulting potential is also known as Yukawa potential in nuclear physics. In the 
derivation of the scheme, we rely on the multipole expansion formula for strictly 
positive A. However, it turns out that the first term of the series representation of 
multipole expansion (with A > 0) corresponds to the case of A = 0. Hence, readers 
can assume that our scheme can be applied to the Laplacian kernel also. 



2 The Multipole Expansion of Multiwavelet Basis Functions 

2.1 The Multipole Expansion 

Denote by (•, •} the L 2 inner product of complex functions on a bounded domain 
in M 3 . Let be a scalar source function supported on the domain. The potential 
generated by cj> (outside of its support) is given by the following multipole series, 

${x) = {G(x,y),ct>{v)) (3) 

oo p 

= £ £ M«k p (\\\x\\)Y*(x) (4) 

p=0 q=—p 

where x = x/\\x\\ is identified to a point on S 2 . The multipole coefficients M$ are 
given by 

M<> = 8\(i P (\\\y\\)Y p q (y),<t>(y)). (5) 
The functions i p and k p are the modified spherical Bessel and Hankel functions, 

i p(0 = y^W/aW ( 6 ) 

h(r) = ^K p+1/2 (r). (7) 

Since i P (r) and k p (r) exhibit exponential growth and decay respectively, an algo- 
rithm based on the above formula is likely to experience an overflow/underflow. 
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To avoid the issue, as suggested in [7] , we replace i p and k p in the formula by their 
scaled forms, 

i$°(Xr)=i(Xr)/X p (8) 
k£°(Xr) = k(Xr)-X p . (9) 

Assuming 0(r) = 1 by appropriate geometric scaling, the appropriate choice of Ao 
is A itself, however, in order to maintain the generality, we keep clear notational 
distinction between them. 

Remark 1 To avoid any confusion, the normalization of spherical harmonics Yp 
used in this paper needs to be clearly stated before we begin any formulation. We 
utilize exactly the same form presented in J7J, that is, 

An obvious advantage of using this form is that Y p ~ q = Yp , hence, in ([5}, Yp 
appears without minus sign in front of q. We can observe, later in this paper, that 
this property provides us with a better symmetry /sparsity pattern of the multipole 
expansion matrix. 



2.2 Symmetries 

The key equation in the above formula is the multipole expansion ([5]). For nota- 
tional simplicity, in this paper, we omit the constant 8 A, which can be multiplied 
afterward. We denote the product of ip° and Yp in (J5j) by Q$; 

Q q p {X,y)=l^{X\\y\\)Y q {y) (11) 

The function Qp should be replaced with the regular solid harmonics when A = 0. 
In the following section, we can observe that the first term of the series represen- 
tation of Qp is the regular solid harmonics. It is obvious that 

Q q (X,ay) = Q q (aX,y). (12) 

The function Qp enjoys two useful symmetries: firstly, from the normalization of 
Yp employed in this paper, it follows that 

Qp q = Qf. (13) 
Secondly, by change of variables <f> — > -k/2 — <f> in Yp , we can obtain 

Qp(X,y 2 ,yi,y3) = i q Qp{X,yi,y2,y?,)- (14) 

As a result, a multipole expansion matrix presented in this paper possesses similar 
symmetries, which we utilize to reduce the number of elements we have to compute. 
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2.3 Multiwavelet Basis Functions 



In this section, we briefly introduce multiwavelet representation of functions. A 
detailed discussion on the subject can be found in [2]. Denote by A; non-negative 
multi-indices and by 4> k multi-dimensional orthonormal polynomials of degree fc, in 
ith dimension. We further assume that the generating functions 4> k are constructed 
by the Cartesian product of 1-d orthonormal polynomials on [—1,1], 



: k 



(y) = J2^M- (15) 



The above 4> k generate the orthonormal multiwavelet basis functions at ar- 
bitrary level n = 0,1,..., and translation characterized by multi-indices I = 
(Zi, . . . , Ij) with Zj = 0, . . . , 2™ — 1 by formula, 

*U*) = l^ din+1) * k{2i2nX - l) - 1) 0nbM (16) 
I elsewhere 

where &(„,{) = Yli=i[ 2 ~ ni i,2~ n {k + 1)]. In this paper, we take [0,l] d (= b by 
definition) as the computational domain. Beware that we assume that the 1-d 
generating functions <f) ki are orthogonal polynomials defined on [—1,1] (not on 
[0, 1]). This choice of unshifted orthogonal polynomials as the generating functions 
is to simplify the notations in multipole related formulae; we have to evaluate 
multipole expansions with respect to the center of each 6( n ,i)- 

Remark 2 The term orthogonal polynomials can be a source of confusion, which 
we need to clarify before we present any related formula. By the term, we mean a 
sequence of polynomials 4> k of degree k orthogonal to each other with respect to an 
underlying weighting function (as in "orthogonal polynomials and quadratures" ) . 
Since a non-trivial weighting function loses its meaning under scaling, readers 
may think 4> k a synonym of (normalized) Legendre polynomial of degree k. This 
limitation of generating functions to orthogonal polynomials greatly simplifies the 
resulting formulae and makes the conversion matrix more sparse. 

Remark 3 There can be different choices of polynomial basis which are mutually 
orthogonal such as the interpolating basis presented in [2] . Conversion between them 
and Legendre-generated basis is not complicated. The advantage (by symmetry 
and sparsity) of using orthogonal polynomials exceeds the additional cost of basis 
conversion. 



2.4 The Multipole Expansion of <j> k n n 

For any p = 0, 1, . . . and q = —p, . . . ,p, define E^' 9 \n, A) by 

Et q) (ri,x) = (Q q p(^v- c (n < i))A h n,i)(y)) hn>l) (17) 
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where ci n! i) is the center of &( n ,i)- The above equation gives the multipole coeffi- 
cient Mp (without 8 A) of (0 with respect to the center, £7 n n. The inner product 
can be scaled and translated to the standard domain [— 1, l] 3 , 

E { k P ' q \n,\) = <<?I(A, 2~ ( - n+1) y),<t> k (y)) (18) 

v2 

= ^4 M) (M/2") (19) 
V2 



where 



4 M) (0,A n ) = ^3 ■),/}• (20) 

V2 



Thus, we are required to evaluate (|20p for arbitrary An which depends on A and the 
level, n. Viewing (p, q) as a row multi-index and k as a column multi-index, E^'^ 
acts as the conversion (multipole expansion) matrix for multiwavelet represented 
functions; let s k n ^ be multiwavelet coefficients for a fixed (n, I). The multipole 

expansion centered at C( n ,i) °f the function s k n « <f) k n n is given by the matrix- 
vector multiplication, 

k 

The same matrix can be used for the conversion from a local expansion to its 
multiwavelet representation. Consider a local expansion with the coefficients Lp, 

oo p 

*0"0 = E E LlQl(\,x). (21) 

p=0 q=—p 

The projection of $ on to the span of the multiwavelet basis is, from the orthonor- 
mality, given by 

= (*, **»,!)> = E ^P("> A ) L p- ( 22 ) 
(p.g) 

i.e., by the multiplication with the conjugate transpose of E^' q \ 



2.5 Symmetries 

Recall the symmetries of Qp. The following two conditions are the immediate 
consequences of (|13p and (|14[) . 

^-")(n,A)=£p)(n,A) (23) 
and 

<i, fe 3)^ A )=H) 9 <i, fc 3)(-' A )- (24) 

In a later section, we will show that, depending on (p,q,k), (I) E^' q ^ is either 
real or pure imaginary and (2) has pre-determined sparsity patterns. Combined 
with the above symmetries, we recommend the following storage for the multipole 
expansion matrices. For each level n, we compute E k p ' q - > for q > and ki > k%, 
and store non-negative q portion of the matrices in two sparse matrices, one for 
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real and the other for imaginary. Separated storage is simple and advantageous 
in the implementation; (i) since each element is either real or imaginary, the two 
sparse matrices have disjoint index sets; it does not require any additional stor- 
age or computation cost due to duplicated indices, (ii) For rows with q < 0, the 
multiplication can be omitted; for a complex vector s* n t y 

k k 

with the negative sign for imaginary matrix. 



2.6 Numerical Issues 

There are three major numerical issues which make the evaluation of E^' q \n, A) 
non-trivial. 

(1) Non-homogeneity of Qp(\, •): Unlike regular solid harmonics, Qp are not 
homogeneous. Since we cannot extract A out of the integral, we have to build 
different E^' q \n,\) depending on A and n, which rules out the possibility of 
utilizing a precomputed table. Since they are not even polynomials, there is 
no simple quadrature which produces the exact integral. Any naive approach 
using adaptive quadrature becomes impractical for the following reasons. 

(2) Rapid growth of Qp(\, •): The function i p (\r) grows exponentially. Scaling 
by using ip° (Ar) helps preventing overflow. However, the function still exhibits 
near singularity for large A. 

(3) Oscillating behavior of <fi h : Although they are polynomials, <j> k have all their 
zeros on (—1, 1). Hence, for large k, an adaptive integrator will encounter with 
highly oscillating integrands. 

From the above characteristics, any adaptive integration requires a large amount 
of computation, or simply it fails to converge especially for highly oscillating cases. 
Beware that, to build N conversion matrices (up to depth level N — 1) for the mul- 
tipole expansion (up to degree P) of functions represented by multiwavelets (of 
degree up to K), we have to compute 0(N x P 2 x K 3 ) elements! 

Our approach begins with rewriting the Qp in a series form. Each term in the 
resulting power series involves a regular solid harmonics weighted by an even power 
of ||a;||, hence, a homogeneous polynomial in R 3 . The Cartesian expansion of this 
polynomial can be explicitly written and its projection on multiwavelet basis can 
be obtained exactly without using any numerical quadrature. 



3 The Series Form 

In this section, we present a series representation of the multipole expansion ma- 
trix, E^' q \n, A). We begin with the identity, 

J^L 1 / r \2m+a 

= £ m ir {m + a + i) (5) ■ ^ 

m=0 
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Utilizing T(m + 1/2) = (2m)!/(4 m ml), we obtain 

/ \ \ p 00 1 / \ \ 2m 

= g G m l( 2 m + 2p + 1)ll (£) ' < 26 > 

Therefore Qp can be written in power series given by 

/ \ \ p 00 i /\ 2 \ m 

QU^)={±) E ml , * + 1V , (t) NI 2m+P W ( 2 7) 



..J (2m + 2p+ 1) 

m=0 v ^ ; 



Define iip, m by 



fl?,m(x) = CyM' 1 \\xf m+p Y«(x) (28) 



where 



Notice that the function Rp^ m is a regular solid harmonics multiplied by ||a;|| 2m , 
hence, a homogeneous polynomial of degree (2m+p). The factor Cy(p, simpli- 
fies the Cartesian expansion of Rp^ m , which we introduced in the following section. 

From the above series representation, iw p (0, A n ) is given by 
v2 

00 / A 2 \ m 

= C E (\n/\o,P,q) J2 Am & T MP,?.*) (30) 



rn=0 



where 



C E (\ n /\ , P , q ) = ( ^Y, (31) 

V2 (2p+l)l\ V-V 



A ( v ) — ( 2 P +1 ) !! (32 ) 
Am(p) " m!(2m + 2p+l)!!' ^ 



I m (p,q,k) = / Rl. m {x)<f) (x)dx. (33) 
•/[-MP 

The factor (2p+ 1)!! in j4 m (p) is added for the normalization, Ao(p) = 1. 

In (|30p . A„ is now taken out of the integral. We will observe that the A- 
independent term, I m {p, q, k), can be further reduced to a finite sum of products of 
1-d integrals with two integer parameters, namely, ij,. We can construct l\, exactly 
without using any numerical quadrature via the recurrence relation of orthogonal 
polynomials. Our strategy is to tabulate P k and use the table to evaluate the series 
(|30|) for various A, n, p, q, and k. 
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3.1 Properties of I m and the Convergence Criterion 

Most of the properties of I m presented in this section will be explained in detail 
in 21 For a more comprehensive presentation, we think it would be more appro- 
priate to discuss the behavior of (|30[) prior to the presentation of detailed formulae 
for I m - Followings are the summary of the relevant properties: 

1. I m (p,q, k) is either real or pure imaginary depending only on k. 

2. Im{p, q,k) = if 2m < k x + k y + k z — p. 

3. Sign of I m (p,q,k) is determined by q only and is independent of m. 

4. \Im(p, q, k)\ is monotonically increasing as m increases. 

5. I rn+ i(p,q,k)/Im(p,q,k) -> 3 as m -> oo. 

Remark 4 Property (5) can be supported by the following estimate: Since \Cy~ x Yp \ 



Thus, the series consists of two parts: A m decreasing factorially and (A 2 /8) m I m 
which behaves asymptotically ~ (3A 2 /8) m . Their product C m = A m (X 2 /8) m I m has 
a fixed sign for a fixed (p, q, k) independently of m. Hence, the partial sum of the 
series increases (decreases) monotonically to the upper (lower) bound which is 
potentially huge in the absolute sense. The non-alternating feature of the series 
suppresses any necessity of considerations of cancellation errors, and suggest the 
following simple convergence criterion: for given absolute and relative tolerances 
t a and e r , stop the summation if 



We can numerically observe that the number of terms to convergence M is 
O(A) in a conservative estimation. For example, for e r = 1CF 16 , M ~ A and slightly 
smaller if A is large; e.g., when A ~ 300, M ~ 200. The condition (2) combined 
with the convergence criterion provides us with additional sparsity of E k p ' q ^; if 
2M < k x + k y + k z — p, the corresponding E^'^ can be considered to be zero. 

Remark 5 We represent E^' q ^ like a special function of A with exponential growth. 
The number of terms M can grow indefinitely as A increases. Although, in many 
practical applications, A are quite limited and \ n = A/2" decreases as the depth of 
the multiwavelet representation increases, a more complete algorithm requires an 
asymptotic expansion of E^'^ with respect to A. Yet, we haven't found a closed 
formula for the asymptotic expansion, which is an on-going work. 



\Pl q \cos6)e iq ^\ < 1, 





(34) 
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4 The Formula for I m and the Sparsity Pattern 



In this section, we present an explicit Cartesian expansion form of Rp, m in I m . 
Each term can be written as a product of 1-d integrals which can be evaluated 
exactly by the recurrence relations of the orthogonal polynomials <f> k . We begin 
with the series form of the spherical harmonics. With the Rodrigues' formula, the 
associated Legendre functions Pp in Yp can be written as 



p v ' 2Pp\ y ' dzP+\i\ K 



(zl^ (1 _, 2) M/ 2 ± (-1H2P-2,)! (35) 



2p y ' ^ i/!(p-z/)!(p-|g|-2i/)! 
Hence, using notations x = (x,y,z), r = \\x\\, and s = sign(g), 

|g| 



Rp, m (*) = r 2m+P P P ql (z/r) 



x — sty 
\fr 2 — z 2 



" 2P {X SlV) A, v\(p-v)!( P -\ q \-2v)! r Z (M) 

By expanding (a; — siy)' 9 ' and r 2< ^ m+v \ we obtain 

= (si) M 
\9\ L^J m +u 

M=0 l/=0 a=0 0=0 



where the coefficients are given by 



{-l) v {2p-2v- 1)!! 
2" v\{p-q-2v)\ 



(38) 
(39) 



= Q (41) 

where we use the definition, (—1)!! = 0!! = 1. 
4.1 The Formula 

From (|37|) . we obtain our final formula for I m : 

I m (p,q,k) = (sz)^ I^(p,q,k) + (si)^ + 1 l£\p,q,k) (42) 
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where 

[llij ^ p-|<?I j 

l£ ) (p,q,k) = X)(-l)"o2^ £ K 

■ J2 °»° Ti p - lql+2m - 2a) J2 d af 3 t^ +2p ~^ iiT' 2 ^, (43) 

a=0 0=0 



i£\ P ,q,k)= (-ir^+i £ ^ 

m-\-v a 
■ J2 C - /If l9l + 2 " l - 2tt) £ d Q0 /tH+2/3-2^-1) J <2a-V+2 M+ l) i {u) 
a=0 0=0 

and 

n=/ c'Ack. (45) 

Note that /m'' and are real functions and, with factors (si)' 9 ' and (si)' 9 ' +1 
respectively, they determine real and imaginary parts of I m separately. The above 
representation of I m by two separate parts Im^ and 1$ is to signify the following 
very useful fact: at least, one of 1$ and iff vanishes for any (p, q, ft) independently 
of m, which implies that 

E k P ' q ^ is either a real or a pure imaginary. 

Moreover, depending on the parameter (p,q,k), many of I m vanish, which results 
in the nice sparsity of the multipole expansion matrix. These properties are the 
immediate consequence of the following properties of orthogonal polynomials. 



(1) Oddity Any orthogonal polynomial <f> with symmetric domain and weight is 
even (odd) if the degree k is even (odd). Hence, 

f k = if (J + ft) is odd. 

Since m always appears in the equation with the factor of 2, any consequence 
of the oddity condition is m-independent; the resulting sparsity of E k P ' q ^ is pre- 
determined by (p, q, k) only (independently of the level n and A). We can observe 
that 

(a) li 2) = if q = 0. 

(b) iff = if k x is odd or (\q\ + k y ) is odd or (p + \q\ + k z ) is odd. 

( 2^ 

(c) I m = if k x is even or (\q\ + k y ) is even or (p + \q\ + k z ) is odd. 
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Therefore, 



Suppose Ej?' q ^ ^= from the above test. Then, the oddity of k x must be the same 
as the oddity of (\q\ + k y ), which results in 



if 



{k z +p + \q\) is odd 
(k x + k y + \q\) is odd 

q = and at least one of k x and k y is odd 



(46) 



r(2) 
m 



if k x is odd 
if k x is even 



(— 1) La J sign(g) i 
(-1)* 



if fey is odd 
if fey is even 



(47) 



Notice, -E^' 9 ^ is real (imaginary) if k y is even (odd). 

The following table illustrates the sparsity of the multipole expansion matrix 
for parameters: < p < 10, < q < p, and < k x , y ,z < 10. We can observe that 
about a quarter of elements are non-zeroes. (See Table [TJ) 



total elements 


real non-zeroes 


imaginary non-zeroes 


87846 


12186 (13.9%) 


8450 (9.6%) 



Table 1 A-independent sparsity estimated by the oddity condition. 



(2) Moment condition Recall the moment conditions satisfied by orthogonal 
polynomials. 

7{ = if I < k. 

Consider a term with a fixed set of indices (p,q,k, n,v,a, f3) in (|42[) . The term 
vanishes if 

k z > p - \q\ + 2m - 2a or k y > \q\ + 2/3 - /i or k x > 2a - 2/3 + fi 
which is true if 

k x + k y + k z > p + 2m. 

Thus, 

l m (p,q,k)=0 if 2m < k x + k y + k z -p (48) 

This condition is m-dependent, hence, cannot be used to pre-determine the sparsity 
pattern of E^' 9 ^ . However, it still can affect the sparsity for a given A; suppose that 
the convergence criterion (|34p is satisfied at M for 2M < k x + k y + k z — p. Then, 
the corresponding E^' q ^ (n, A) is effectively zero. Since M ~ A and E^' q \n, A) = 

constant ■ E^ ,q) {Q, A/2"), the multipole expansion matrix becomes more sparse as 
A decreases and as n increases. Table [2] shows the number of vanishing elements 
(for n = 0) among those predicted to be non-zeroes by the oddity condition. With 
parameters, < p < 10, < q < p, and < k x , y , z < 10, the number of elements is 
87846 and the numbers of non-zero real and imaginary elements (predicted by the 
oddity condition) are 12186 and 8450 respectively (same as the above example). 
Tolerances are e a = e r = 10~ 16 . 
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A 


additional real zeroes 


additional imaginary zeroes 


1 


9567 


6679 


2 


8813 


6154 


4 


7478 


5235 


6 


6340 


4439 


8 


5439 


3775 


10 


4630 


3203 


50 


1309 


851 


100 


1301 


848 


200 


1273 


835 


300 


1251 


824 



Table 2 A-dcpcndcnt sparsity estimated from the moment condition. 



The additional sparsity decreases as A, hence M, increases. There are two 
factors which controls M and, hence, the the additional sparsity - the absolute 
tolerance e a and the relative tolerance e r . Among them, the contribution of e r 
decreases rapidly and becomes quite negligible when A 3> max fcj. However, the 
contribution of e a is persistent. From the table, we can observe that the additional 
sparsities by e a are ~ 1200 for the real matrix and ~ 800 for the imaginary matrix. 

Also, the moment condition enhances the computational efficiency slightly. 
We can simply skip the evaluation of I m if m < \{k x + k y + k z — p)/2~|. Like 
the sparsity by e r , the effect decreases quite rapidly as A increases. However, in 
practical applications of fast multipole and multiwavelet representation, the levels 
of terminal boxes where we need to perform the expansion is likely to be high. 
Hence, the additional sparsity by e r should not be considered insignificant. 

4.2 The Laplacian Kernel (A = 0) 

In this case, Qp becomes simply the regular solid harmonics. The corresponding 
E^' qS> (n, 0) is just the first term (m = 0) of the series form (|30|l with an appropriate 
adjustment of the constant factor. In this case, the moment condition becomes 

E k p ' 9) =0 if p<k x + k y + k z 

which results in a more sparse multipole expansion matrices. With the same con- 
dition as the previous examples, maxp = maxfc.; = 10, the number of non-zero 
elements are given in Table [3] 



total elements real non-zeroes imaginary non-zeroes 
87846 1512 (1.72%) 1001 (1.14%) 



Table 3 Sparsity of the case with A = 0. 



We can observe that the resulting matrices are very sparse - only less than 
3% of total elements are non-zeroes. This example illustrates the efficiency of the 
multipole expansion on multiwavelet representations based on orthogonal (almost 
synonymously in this paper, Legendre) polynomials. 
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4.3 Recurrence Relations for ij. 

Finally, we present the algorithm to build the table of l\. required for the evaluation 
of 7 m - Let p < pmax, fcj < fc max , and M < M max (for a given A). Then, the required 
size of table is (2M max +p m ax) X fc max , and the half of the elements are zero by 
the oddity condition. 

Each element if, can be calculated by the identical recurrence relations to those 
of the orthogonal polynomials <f> k . Recall any sequence orthogonal polynomials 
satisfy a three term recurrence relation of the form, 

Z+^KC + ft)/-^^ 1 . (49) 

It immediately follows that 

I k+ 1 = a k T k +1 + p k T k - 7fc 7 fe _!. (50) 

From the oddity condition, if. , 1 ^ if and only if 1^ = 0, and r k = for I < k. 
Hence, for k < I, 



ji+i _ I if (I + k) is even 

' 1 ^A+i+Vk^klk-i if G + is odd 



and 

The recurrence relation can be evaluated from the initial data, 



n = or k ^J k k -_\ (v?;! = o). (52) 



% = ao 1 || ^ and % = a l 1 ^ (53) 

where 0° = ao and c^ 1 = ai£. 

For normalized Legendre polynomials, the coefficients are given by 



k + 1 2k + l , -i /2fe+l 

2fcTTV2fcT3 and ^ 7fe = 2fcTIV2fc^T- (54) 



5 Results and Conclusions 

Most of the implementation is done very faithfully with the formulae presented in 
this paper. The only special treatment is that, in order to suppress the accumula- 
tion of round-off errors affecting the result, we used higher precision floating point 
arithmetics for internal calculations including the table for ij,; for example, to 
generate matrices with 64bit double precision, we employed 80bit long(extended)- 
double arithmetics. The computational cost is governed by the number of terms to 
be added, m, and is not significantly affected by the augmented internal precision. 
By comparing with values obtained by applying adaptive numerical integrator to 
(|15p . we could validate the presented formula. When A or k is only moderately 
large, an adaptive integrator typically fails to converge since the integrand of (U5[) 
becomes near-singular or highly oscillating. Thus, the presented formula can be 
viewed as a reliable way to evaluate (|15|) (or a similar form of integral) when a typ- 
ical numerical quadrature is not applicable due to the near-singularity and/or the 



14 



Jac-Scok Huh 



oscillation of the integrand. It is also observed that the computing time of building 
l\. table is negligible compared to the computing time of E^' q K We summarize 
the contributions of this paper as follows. 

1. We presented a method to build the multipole expansion matrices for functions 
represented by multiwavelets. 

2. The presented method does not involve any numerical quadrature and based 
entirely on a series representation like a special function of A. 

3. The proposed scheme generates highly accurate multipole conversion matrices 
stably and reliably for a wide range of parameters (p, q, k) and A. 
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